/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __LOGISTIC_LOSS_HH__
#define __LOGISTIC_LOSS_HH__

#include <algorithm>
#include <boost/optional.hpp>
#include <cmath>
#include <stdint.h>
#include "loss.hh"

#include "wexpxmx.hh"

namespace flassol
{
  class LogisticLoss : public Loss
    {
      private:
        static LossFactory factory;

        uint64_t bad_labels;
        uint64_t npos;
        uint64_t nneg;

      public:
        explicit
        LogisticLoss () : bad_labels (0), npos (0), nneg (0) { }

        virtual float 
        loss (float p,
              float z)
          {
            if (z != 1 && z != 0)
              {
                ++bad_labels;
              }

            float y = (z > 0) ? 1 : -1;

            if (y == 1) 
              { 
                ++npos; 
              } 
            else
              {
                ++nneg;
              }

            if (y * p > 10)
              {
                // Log[1 + Exp[-x]] = Exp[-x] + O (Exp[-2x])
                
                return expf (-y * p);
              }
            else if (y * p < -10)
              {
                // Log[1 + Exp[-x]] = Log[Exp[-x]] + Log[Exp[x] + 1]
                //                 = -x + O (Exp[x])
                
                return -y * p;
              }
            else
              {
                return logf (1.0 + expf (-y * p));
              }
          }

        virtual float 
        dldp (float p,
              float z)
          {
            float y = (z > 0) ? 1 : -1;

            if (y * p > 13)
              {
                return -y * expf (-y * p);
              }
            else
              {
                return -y / (1.0 + expf (y * p));
              }
          }

        virtual float 
        invariant_update (float p,
                          float wdotx,
                          float adotb,
                          float z,
                          float eta,
                          float anormplusbnorm,
                          float xnorm,
                          float onedotone)
          {
            // TODO: handle dyadic case
            (void) wdotx;
            (void) adotb;
            (void) anormplusbnorm;
            (void) onedotone;

            float y = (z > 0) ? 1 : -1;

            // The actual learning rate is (eta / norm)
            //
            // (W (exp (eta + y p + exp (y p))) - eta - exp (y p)) / (y norm)
            // = (G (eta + y p + exp (y p)) + y p) / (y norm)
            // = (y G (eta + y p + exp (y p)) + p) / norm
            //
            // where G (x) = W (exp (x)) - x

            return (y * wexpxmx (eta + y * p + expf (y * p)) + p) / xnorm;
          }

        virtual std::ostream& 
        report (std::ostream& o,
                LogLevel      level)
          {
            if (bad_labels > 0 && level <= WARN)
              {
                o << "WARNING: encountered " << bad_labels 
                  << " labels that were not 1 or 0"
                  << std::endl;
                o << "WARNING: positives were mapped to 1"
                  << " non-positives were mapped to 0"
                  << std::endl;
              }

            if (npos + nneg && level <= INFO)
              {
                double p =   static_cast<double> (npos) 
                           / static_cast<double> (npos + nneg);
                double Hp = (npos) ? npos * log (p) : 0;

                Hp += (nneg) ? nneg * log (1.0 - p) : 0;

                o << "unconditional probability of class 1 = " << p << std::endl;
                o << "best constant = ";

                if (p == 0)
                  {
                    o << "-infinity";
                  }
                else if (p == 1)
                  {
                    o << "+infinity";
                  }
                else
                  {
                    o << log (p / (1.0 - p));
                  }

                o << std::endl;

                o << "best constant loss = " << -Hp / (npos + nneg) << std::endl;
              }

            return o;
          }
    };
}

#endif // __LOGISTIC_LOSS_HH__
